Foreword

This presentation borrows intends to be an introduction to the practical aspects of Machine Learning for students of the Computational Methods Uniandes course. The content provided below borrows heavily from Coursera’s Data Science Specialization by Jeff Leek, Roger Peng and Brian Caffo (Johns Hopkins University).

Mathematical Foundations of Prediction

Prediction

• process for prediction = population –> probability and sampling to pick set of data –> split into training and test set –> build prediction function –> predict for new data –> evaluate
• Note: choosing the right dataset and knowing what the specific question is are paramount to the success of the prediction algorithm (GoogleFlu failed to predict accurately when people’s search habits changed)

• components of predictor = question –> input data –> features (extracting variables/characteristics) –> algorithm –> parameters (estimate) –> evaluation
• relative order of importance = question (concrete/specific) > data (relevant) > features (properly extract) > algorithms
• data selection
  • Note: “garbage in = garbage out” –> having the correct/relevant data will decide whether the model is successful
  • data for what you are trying to predict is most helpful
  • more data –> better models (usually)
• feature selection
  • good features –> lead to data compression, retain relevant information, created based on expert domain knowledge
  • common mistakes –> automated feature selection (can yield good results but likely to behave inconsistently with slightly different data), not understanding/dealing with skewed data/outliers, throwing away information unnecessarily
• algorithm selection
  • matter less than one would expect
  • getting a sensible approach/algorithm will be the basis for a successful prediction
  • more complex algorithms can yield incremental improvements
  • ideally interpretable (simple to explain), accurate, scalable/fast (may leverage parallel computation)
• prediction is effectively about trade-offs
  • find the correct balance between interpretability vs accuracy vs speed vs simplicity vs scalability
  • interpretability is especially important in conveying how features are used to predict outcome
  • scalability is important because for an algorithm to be of practical use, it needs to be implementable on large datasets without incurring large costs (computational complexity/time)

In Sample vs Out of Sample Errors

• in sample error = error resulted from applying your prediction algorithm to the dataset you built it with
  • also known as resubstitution error
  • often optimistic (less than on a new sample) as the model may be tuned to error of the sample
• out of sample error = error resulted from applying your prediction algorithm to a new data set
  • also known as generalization error
  • out of sample error most important as it better evaluates how the model should perform
• in sample error < out of sample error
  • reason is over-fitting: model too adapted/optimized for the initial dataset
    • data have two parts: signal vs noise
    • goal of predictor (should be simple/robust) = find signal
    • it is possible to design an accurate in-sample predictor, but it captures both signal and noise
    • predictor won’t perform as well on new sample
  • often times it is better to give up a little accuracy for more robustness when predicting on new data
• example

# load data
library(kernlab); data(spam); set.seed(333)
# picking a small subset (10 values) from spam data set
smallSpam <- spam[sample(dim(spam)[1],size=10),]
# label spam = 2 and ham = 1
spamLabel <- (smallSpam$type=="spam")*1 + 1
# plot the capitalAve values for the dataset with colors differentiated by spam/ham (2 vs 1)
plot(smallSpam$capitalAve,col=spamLabel)

# first rule (over-fitting to capture all variation)
rule1 <- function(x){
    prediction <- rep(NA,length(x))
    prediction[x > 2.7] <- "spam"
    prediction[x < 2.40] <- "nonspam"
    prediction[(x >= 2.40 & x <= 2.45)] <- "spam"
    prediction[(x > 2.45 & x <= 2.70)] <- "nonspam"
    return(prediction)
}
# tabulate results of prediction algorithm 1 (in sample error --> no error in this case)
table(rule1(smallSpam$capitalAve),smallSpam$type)

##          
##           nonspam spam
##   nonspam       5    0
##   spam          0    5

# second rule (simple, setting a threshold)
rule2 <- function(x){
    prediction <- rep(NA,length(x))
    prediction[x > 2.8] <- "spam"
    prediction[x <= 2.8] <- "nonspam"
    return(prediction)
}
# tabulate results of prediction algorithm 2(in sample error --> 10% in this case)
table(rule2(smallSpam$capitalAve),smallSpam$type)

##          
##           nonspam spam
##   nonspam       5    1
##   spam          0    4

# tabulate out of sample error for algorithm 1
table(rule1(spam$capitalAve),spam$type)

##          
##           nonspam spam
##   nonspam    2141  588
##   spam        647 1225

# tabulate out of sample error for algorithm 2
table(rule2(spam$capitalAve),spam$type)

##          
##           nonspam spam
##   nonspam    2224  642
##   spam        564 1171

# accuracy and total correct for algorithm 1 and 2
rbind("Rule 1" = c(Accuracy = mean(rule1(spam$capitalAve)==spam$type),
    "Total Correct" = sum(rule1(spam$capitalAve)==spam$type)),
    "Rule 2" = c(Accuracy = mean(rule2(spam$capitalAve)==spam$type),
    "Total Correct" = sum(rule2(spam$capitalAve)==spam$type)))

##         Accuracy Total Correct
## Rule 1 0.7315801          3366
## Rule 2 0.7378831          3395

Prediction Study Design

• procedures
  1. define error rate (type I/type II)
  2. split data into:
     • training, testing, validation (optional)
  3. pick features from the training set
     • use cross-validation
  4. pick prediction function (model) on the training set
     • use cross-validation
  5. if no validation set
     • apply 1 time to test set
  6. if there is a validation set
     • apply to test set and refine
     • apply 1 time to validation
  • Note: it’s important to hold out an untouched sample to accurately estimate the out of sample error rate
• benchmarks (i.e. set all variables = 0) can help pinpoint/test the model to see what is wrong with the model
• avoid small sample sizes
  • consider binary outcomes (i.e. coin flip)
  • for n = 1, the probability of perfect classification (100% accuracy) is 50%
  • for n = 10, the probability of perfect classification (100% accuracy) is 0.1%
  • so it’s important to have bigger samples so that when you do get a high accuracy, it may actually be a significant result and not just by chance
• example: Netflix rating prediction competition
  • split data between training and “held-out”
    • held-out included probe, quiz and test sets
    • probe is used to test the predictor built from the training dataset
    • quiz is used to realistically evaluate out of sample error rates
    • test is used to finally evaluate the validity of algorithm
  • important to not tune model to quiz set specifically

Sample Division Guidelines for Prediction Study Design

• for large sample sizes
  • 60% training
  • 20% test
  • 20% validation
• for medium sample sizes
  • 60% training
  • 40% test
  • no validation set to refine model (to ensure test set is of sufficient size)
• for small sample sizes
  • carefully consider if there are enough sample to build a prediction algorithm
  • no test/validation sets
  • perform cross validation
  • report caveat of small sample size and highlight the fact that the prediction algorithm has never been tested for out of sample error
• there should always be a test/validation set that is held away and should NOT be looked at when building model
  • when complete, apply the model to the held-out set only one time
• randomly sample training and test sets
  • for data collected over time, build training set in chunks of times
• datasets must reflect structure of problem
  • if prediction evolves with time, split train/test sets in time chunks (known as backtesting in finance)
• subsets of data should reflect as much diversity as possible

Picking the Right Data

• use like data to predict like
• to predict a variable/process X, use the data that’s as closely related to X as possible
• weighting the data/variables by understanding and intuition can help to improve accuracy of prediction
• data properties matter –> knowing how the data connects to what you are trying to measure
• predicting on unrelated data is the most common mistake
  • if unrelated data must be used, be careful about interpreting the model as to why it works/doesn’t work

Types of Errors

• when discussing the outcome decided on by the algorithm, Positive = identified and negative = rejected
  • True positive = correctly identified (predicted true when true)
  • False positive = incorrectly identified (predicted true when false)
  • True negative = correctly rejected (predicted false when false)
  • False negative = incorrectly rejected (predicted false when true)
• example: medical testing
  • True positive = Sick people correctly diagnosed as sick
  • False positive = Healthy people incorrectly identified as sick
  • True negative = Healthy people correctly identified as healthy
  • False negative = Sick people incorrectly identified as healthy

Notable Measurements for Error – Binary Variables

• accuracy = weights false positives/negatives equally
• concordance = for multi-class cases, \[\kappa = \frac{accuracy - P(e)}{1 - P(e)}\] where \[P(e) = \frac{TP+FP}{total} \times \frac{TP+FN}{total} + \frac{TN+FN}{total} \times \frac{FP+TN}{total}\]

• example
  • given that a disease has 0.1% prevalence in the population, we want to know what’s probability of a person having the disease given the test result is positive? the test kit for the disease is 99% sensitive (most positives = disease) and 99% specific (most negatives = no disease)
  • what about 10% prevalence?

Notable Measurements for Error – Continuous Variables

• Mean squared error (MSE) = \[\frac{1}{n} \sum_{i=1}^n (Prediction_i - Truth_i)^2\]
• Root mean squared error (RMSE) - \[\sqrt{\frac{1}{n} \sum_{i=1}^n(Prediction_i - Truth_i)^2}\]
  • in the same units as variable
  • most commonly used error measure for continuous data
  • is not an effective measures when there are outliers
    • one large value may significantly raise the RMSE
• median absolute deviation = \[median(|Prediction_i - Truth_i|)\]

Receiver Operating Characteristic Curves

• predictions for binary classification often are quantitative (i.e. probability, scale of 1 to 10)
  • different cutoffs/threshold of classification (> 0.8 –> one outcome) yield different results/predictions
  • Receiver Operating Characteristic curves are generated to compare the different outcomes

• ROC Curves
  • x-axis = 1 - specificity (or, probability of false positive)
  • y-axis = sensitivity (or, probability of true positive)
  • points plotted = cutoff/combination
  • areas under curve = quantifies whether the prediction model is viable or not
    • higher area –> better predictor
    • area = 0.5 –> effectively random guessing (diagonal line in the ROC curve)
    • area = 1 –> perfect classifier
    • area = 0.8 –> considered good for a prediction algorithm

• example
  • each point on the graph corresponds with a specificity and sensitivity

Cross Validation

• procedures
  1. split training set into sub-training/test sets
  2. build model on sub-training set
  3. evaluate on sub-test set
  4. repeat and average estimated errors
• result
  • we are able to fit/test various different models with different variables included to the find the best one on the cross-validated test sets
  • we are able to test out different types of prediction algorithms to use and pick the best performing one
  • we are able to choose the parameters in prediction function and estimate their values
  • Note: original test set completely untouched, so when final prediction algorithm is applied, the result will be an unbiased measurement of the out of sample accuracy of the model
• approaches
  • random subsampling
  • K-fold
  • leave one out
• considerations
  • for time series data data must be used in “chunks”
    • one time period might depending all time periods previously (should not take random samples)
  • if you cross-validate to pick predictors, the out of sample error rate may not be the most accurate and thus the errors should still be measured on independent data

Random Subsampling

* a randomly sampled test set is subsetted out from the original training set * the predictor is built on the remaining training data and applied to the test set * the above are three random subsamplings from the same training set * considerations - must be done without replacement - random sampling with replacement = bootstrap + underestimates of the error + can be corrected, but it is complicated (0.632 Bootstrap)

K-Fold

* break training set into \(K\) subsets (above is a 3-fold cross validation) * build the model/predictor on the remaining training data in each subset and applied to the test subset * rebuild the data \(K\) times with the training and test subsets and average the findings * considerations
- larger k = less bias, more variance - smaller k = more bias, less variance

Leave One Out

* leave out exactly one sample and build predictor on the rest of training data * predict value for the left out sample * repeat for each sample

caret Package (tutorial)

• core functionality
  • preprocessing/cleaning data –> preProcess()
  • cross validation/data splitting –> createDataPartition(), createResample(), createTimeSlices()
  • train algorithms on training data and apply to test sets –> train(), predict()
  • model comparison (evaluate the accuracy of model on new data) –> confusionMatrix()
• machine learning algorithms in caret package
  • linear discriminant analysis
  • regression
  • naive Bayes
  • support vector machines
  • classification and regression trees
  • random forests
  • boosting
  • many others
• caret provides uniform framework to build/predict using different models
  • create objects of different classes for different algorithms, and caret package allows algorithms to be run the same way through predict() function

Data Slicing

• createDataPartition(y=data$var, times=1, p=0.75, list=FALSE) –> creates data partitions using given variable
  • y=data$var = specifies what outcome/variable to split the data on
  • times=1 = specifies number of partitions to create (number of data splitting performed)
  • p=0.75 = percent of data that will be for training the model
  • list=FALSE = returns a matrix of indices corresponding to p% of the data (training set)
    • Note: matrix is easier to subset the data with, so list = FALSE is generally what is used
    • list=TRUE = returns a list of indices corresponding to p% of the data (training set)
  • the function effectively returns a list of indexes of the training set which can then be leveraged to subset the data
    • training<-data[inTrain, ] = subsets the data to training set only
    • testing<-data[-inTrain, ] = the rest of the data set can then be stored as the test set
  • example

# load packages and data
library(caret)
# create training set indexes with 75% of data
inTrain <- createDataPartition(y=spam$type,p=0.75, list=FALSE)
# subset spam data to training
training <- spam[inTrain,]
# subset spam data (the rest) to test
testing <- spam[-inTrain,]
# dimension of original and training dataset
rbind("original dataset" = dim(spam),"training set" = dim(training))

##                  [,1] [,2]
## original dataset 4601   58
## training set     3451   58

• createFolds(y=data$var, k=10, list=TRUE, returnTrain=TRUE) = slices the data in to \(k\) folds for cross validation and returns \(k\) lists of indices
  • y=data$var = specifies what outcome/variable to split the data on
  • k=10 = specifies number of folds to create (See K Fold Cross Validation)
    • each training set has approximately has \(\frac{k-1}{k} \%\) of the data (in this case 90%)
    • each training set has approximately has \(\frac{1}{k} \%\) of the data (in this case 10%)
  • list=TRUE = returns k list of indices that corresponds to the cross-validated sets
    • Note: the returned list conveniently splits the data into k datasets/vectors of indices, so list=TRUE is generally what is used
    • when the returned object is a list (called folds in the case), you can use folds[[1]][1:10] to access different elements from that list
    • list=FALSE = returns a vector indicating which of the k folds each data point belongs to (i.e. 1 - 10 is assigned for each of the data points in this case)
      • Note: these group values corresponds to test sets for each cross validation, which means everything else besides the marked points should be used for training
  • [only works when list=T] returnTrain=TRUE = returns the indices of the training sets
    • [default value when unspecified]returnTrain=FALSE = returns indices of the test sets
  • example

# create 10 folds for cross validation and return the training set indices
folds <- createFolds(y=spam$type,k=10,list=TRUE,returnTrain=TRUE)
# structure of the training set indices
str(folds)

## List of 10
##  $ Fold01: int [1:4141] 1 2 3 4 6 7 8 9 10 12 ...
##  $ Fold02: int [1:4140] 1 2 3 4 5 7 8 9 10 11 ...
##  $ Fold03: int [1:4141] 1 3 4 5 6 7 8 9 10 11 ...
##  $ Fold04: int [1:4141] 2 3 4 5 6 7 8 9 10 11 ...
##  $ Fold05: int [1:4141] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Fold06: int [1:4141] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Fold07: int [1:4141] 1 2 3 4 5 6 8 9 11 12 ...
##  $ Fold08: int [1:4141] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Fold09: int [1:4141] 1 2 4 5 6 7 10 11 12 13 ...
##  $ Fold10: int [1:4141] 1 2 3 5 6 7 8 9 10 11 ...

# return the test set indices instead
# note: returnTrain = FALSE is unnecessary as it is the default behavior
folds.test <- createFolds(y=spam$type,k=10,list=TRUE,returnTrain=FALSE)
str(folds.test)

## List of 10
##  $ Fold01: int [1:460] 15 16 18 40 45 62 68 81 82 102 ...
##  $ Fold02: int [1:459] 1 41 55 58 67 75 117 123 151 175 ...
##  $ Fold03: int [1:461] 3 14 66 69 70 80 90 112 115 135 ...
##  $ Fold04: int [1:460] 5 19 25 65 71 83 85 88 91 93 ...
##  $ Fold05: int [1:460] 6 10 17 21 26 56 57 104 107 116 ...
##  $ Fold06: int [1:459] 7 8 13 39 52 54 76 89 99 106 ...
##  $ Fold07: int [1:461] 4 23 27 29 32 33 34 38 49 51 ...
##  $ Fold08: int [1:460] 2 9 30 31 36 37 43 46 47 48 ...
##  $ Fold09: int [1:461] 12 20 24 44 53 59 60 64 84 98 ...
##  $ Fold10: int [1:460] 11 22 28 35 42 61 72 86 92 118 ...

# return first 10 elements of the first training set
folds[[1]][1:10]

##  [1]  1  2  3  4  6  7  8  9 10 12

• createResample(y=data$var, times=10, list=TRUE) = create 10 resamplings from the given data with replacement
  • list=TRUE = returns list of n vectors that contain indices of the sample
    • Note: each of the vectors is of length of the data, and contains indices
  • times=10 = number of samples to create

# create 10 resamples
resamples <- createResample(y=spam$type,times=10,list=TRUE)
# structure of the resamples (note some samples are repeated)
str(resamples)

## List of 10
##  $ Resample01: int [1:4601] 1 4 4 4 7 8 12 13 13 14 ...
##  $ Resample02: int [1:4601] 3 3 5 7 10 12 12 13 13 14 ...
##  $ Resample03: int [1:4601] 1 2 2 3 4 5 8 10 11 12 ...
##  $ Resample04: int [1:4601] 1 3 3 4 7 8 8 9 10 14 ...
##  $ Resample05: int [1:4601] 2 4 5 6 7 7 8 8 9 12 ...
##  $ Resample06: int [1:4601] 3 6 6 7 8 9 12 13 13 14 ...
##  $ Resample07: int [1:4601] 1 2 2 5 5 6 7 8 9 10 ...
##  $ Resample08: int [1:4601] 2 2 3 4 4 7 7 8 8 9 ...
##  $ Resample09: int [1:4601] 1 4 7 8 8 9 12 13 15 15 ...
##  $ Resample10: int [1:4601] 1 3 4 4 7 7 9 9 10 11 ...

• createTimeSlices(y=data, initialWindow=20, horizon=10) = creates training sets with specified window length and the corresponding test sets
  • initialWindow=20 = number of consecutive values in each time slice/training set (i.e. values 1 - 20)
  • horizon=10 = number of consecutive values in each predict/test set (i.e. values 21 - 30)
  • fixedWindow=FALSE = training sets always start at the first observation
    • this means that the first training set would be 1 - 20, the second will be 1 - 21, third 1 - 22, etc.
    • but the test sets are still like before (21 - 30, 22 - 31, etc.)

# create time series data
tme <- 1:1000
# create time slices
folds <- createTimeSlices(y=tme,initialWindow=20,horizon=10)
# name of lists
names(folds)

## [1] "train" "test"

# first training set
folds$train[[1]]

##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

# first test set
folds$test[[1]]

##  [1] 21 22 23 24 25 26 27 28 29 30

Training Options (tutorial)

• train(y ~ x, data=df, method="glm") = function to apply the machine learning algorithm to construct model from training data

# returns the arguments of the default train function
args(train.default)

## function (x, y, method = "rf", preProcess = NULL, ..., weights = NULL, 
##     metric = ifelse(is.factor(y), "Accuracy", "RMSE"), maximize = ifelse(metric == 
##         "RMSE", FALSE, TRUE), trControl = trainControl(), tuneGrid = NULL, 
##     tuneLength = 3) 
## NULL

• train function has a large set of parameters, below are the default options
  • method="rf" = default algorithm is random forest for training a given data set; caret contains a large number of algorithms
    • names(getModelInfo()) = returns all the options for method argument
    • list of models and their information can be found here
  • preProcess=NULL = set preprocess options (see Preprocessing)
  • weights=NULL = can be used to add weights to observations, useful for unbalanced distribution (a lot more of one type than another)
  • metric=ifelse(is.factor(y), "Accuracy", "RMSE") = default metric for algorithm is Accuracy for factor variables, and RMSE, or root mean squared error, for continuous variables
    • Kappa = measure of concordance (see Notable Measurements for Error – Binary Variables)
    • RSquared can also be used here as a metric, which represents R² from regression models (only useful for linear models)
  • maximize=ifelse(metric=="RMSE", FALSE, TRUE) = the algorithm should maximize accuracy and minimize RMSE
  • trControl=trainControl() = training controls for the model, more details below
  • tuneGrid=NULL
  • tuneLength=3

# returns the default arguments for the trainControl object
args(trainControl)

## function (method = "boot", number = ifelse(grepl("cv", method), 
##     10, 25), repeats = ifelse(grepl("cv", method), 1, number), 
##     p = 0.75, initialWindow = NULL, horizon = 1, fixedWindow = TRUE, 
##     verboseIter = FALSE, returnData = TRUE, returnResamp = "final", 
##     savePredictions = FALSE, classProbs = FALSE, summaryFunction = defaultSummary, 
##     selectionFunction = "best", preProcOptions = list(thresh = 0.95, 
##         ICAcomp = 3, k = 5), index = NULL, indexOut = NULL, timingSamps = 0, 
##     predictionBounds = rep(FALSE, 2), seeds = NA, adaptive = list(min = 5, 
##         alpha = 0.05, method = "gls", complete = TRUE), allowParallel = TRUE) 
## NULL

• trainControl creates an object that sets many options for how the model will be applied to the training data
  • Note: the default values are listed below but you can use them to set the parameters to your discretion
  • method="boot" =
    • "boot" = bootstrapping (drawing with replacement)
    • "boot632" = bootstrapping with adjustment
    • "cv" = cross validation
    • "repeatedcv" = repeated cross validation
    • "LOOCV" = leave one out cross validation
  • number=ifelse(grepl("cv", method),10, 25) = number of subsamples to take
    • number=10 = default for any kind of cross validation
    • number=25 = default for bootstrapping
    • Note: number should be increased when fine-tuning model with large number of parameter
  • repeats=ifelse(grepl("cv", method), 1, number) = numbers of times to repeat the subsampling
    • repeats=1 = default for any cross validation method
    • repeats=25 = default for bootstrapping
  • p=0.75 = default percentage of data to create training sets
  • initialWindow=NULL, horizon=1, fixedWindow=TRUE = parameters for time series data
  • verboseIter=FALSE = print the training logs
  • returnData=TRUE, returnResamp = “final”,
  • savePredictions=FALSE = save the predictions for each resample
  • classProbs=FALSE = return classification probabilities along with the predictions
  • summaryFunction=defaultSummary = default summary of the model,
  • preProcOptions=list(thresh = 0.95, ICAcomp = 3, k = 5) = specifies preprocessing options for the model
  • predictionBounds=rep(FALSE, 2) = specify the range of the predicted value
    • for numeric predictions, predictionBounds=c(10, NA) would mean that any value lower than 10 would be treated as 10 and no upper bounds
  • seeds=NA = set the seed for the operation
    • Note: setting this is important when you want to reproduce the same results when the train function is run
  • allowParallel=TRUE = sets for parallel processing/computations

Plotting Predictors (tutorial)

• it is important to only plot the data in the training set
  • using the test data may lead to over-fitting (model should not be adjusted to test set)
• goal of producing these exploratory plots = look for potential outliers, skewness, imbalances in outcome/predictors, and explainable groups of points/patterns
• featurePlot(x=predictors, y=outcome, plot="pairs") = short cut to plot the relationships between the predictors and outcomes

# load relevant libraries
library(ISLR); library(ggplot2);
# load wage data
data(Wage)
# create training and test sets
inTrain <- createDataPartition(y=Wage$wage,p=0.7, list=FALSE)
training <- Wage[inTrain,]
testing <- Wage[-inTrain,]
# plot relationships between the predictors and outcome
featurePlot(x=training[,c("age","education","jobclass")], y = training$wage,plot="pairs")

• qplot(age, wage, color=eduction, data=training) = can be used to separate data by a factor variable (by coloring the points differently)
  • geom_smooth(method = "lm") = adds a regression line to the plots
  • geom=c("boxplot", "jitter") = specifies what kind of plot to produce, in this case both the boxplot and the point cloud

# qplot plus linear regression lines
qplot(age,wage,colour=education,data=training)+geom_smooth(method='lm',formula=y~x)

• cut2(variable, g=3) = creates a new factor variable by cutting the specified variable into n groups (3 in this case) based on percentiles
  • Note: cut2 function is part of the Hmisc package, so library(Hmisc) must be run first
  • this variable can then be used to tabulate/plot the data
• grid.arrange(p1, p2, ncol=2) = ggplot2 function the print multiple graphs on the same plot
  • Note: grid.arrange function is part of the gridExtra package, so library(gridExtra) must be run first

# load Hmisc and gridExtra packages
library(Hmisc);library(gridExtra);
# cute the wage variable
cutWage <- cut2(training$wage,g=3)
# plot the boxplot
p1 <- qplot(cutWage,age, data=training,fill=cutWage,
      geom=c("boxplot"))
# plot boxplot and point clusters
p2 <- qplot(cutWage,age, data=training,fill=cutWage,
      geom=c("boxplot","jitter"))
# plot the two graphs side by side
grid.arrange(p1,p2,ncol=2)

• table(cutVariable, data$var2) = tabulates the cut factor variable vs another variable in the dataset
• prop.table(table, margin=1) = converts a table to a proportion table
  • margin=1 = calculate the proportions based on the rows
  • margin=2 = calculate the proportions based on the columns

# tabulate the cutWage and jobclass variables
t <- table(cutWage,training$jobclass)
# print table
t

##                
## cutWage         1. Industrial 2. Information
##   [ 20.1, 92.7)           445            256
##   [ 92.7,119.1)           376            347
##   [119.1,318.3]           269            409

# convert to proportion table based on the rows
prop.table(t,1)

##                
## cutWage         1. Industrial 2. Information
##   [ 20.1, 92.7)     0.6348074      0.3651926
##   [ 92.7,119.1)     0.5200553      0.4799447
##   [119.1,318.3]     0.3967552      0.6032448

• qplot(var1, color=var2, data=training, geom="density") = produces density plot for the given numeric and factor variables
  • effectively smoothed out histograms
  • provides for easy overlaying of groups of data
    • break different variables up by group and see how outcomes change between groups

# produce density plot
qplot(wage,colour=education,data=training,geom="density")

Preprocessing (tutorial)

• some predictors may have strange distributions (i.e. skewed) and may need to be transformed to be more useful for prediction algorithm
  • particularly true for model based algorithms –> naive Bayes, linear discriminate analysis, linear regression
• centering = subtracting the observations of a particular variable by its mean
• scaling = dividing the observations of a particular variable by its standard deviation
• normalizing = centering and scaling the variable –> effectively converting each observation to the number of standard deviations away from the mean
  • the distribution of the normalized variable will have a mean of 0 and standard deviation of 1
  • Note: normalizing data can help remove bias and high variability, but may not be applicable in all cases
• Note: if a predictor/variable is standardized when training the model, the same transformations must be performed on the test set with the mean and standard deviation of the train variables
  • this means that the mean and standard deviation of the normalized test variable will NOT be 0 and 1, respectively, but will be close
  • transformations must likely be imperfect but test/train sets must be processed the same way
• train(y~x, data=training, preProcess=c("center", "scale")) = preprocessing can be directly specified in the train function
  • preProcess=c("center", "scale") = normalize all predictors before constructing model
• preProcess(trainingData, method=c("center", "scale") = function in the caret to standardize data
  • you can store the result of the preProcess function as an object and apply it to the train and test sets using the predict function

# load spam data
data(spam)
# create train and test sets
inTrain <- createDataPartition(y=spam$type,p=0.75, list=FALSE)
training <- spam[inTrain,]
testing <- spam[-inTrain,]
# create preProcess object for all predictors ("-58" because 58th = outcome)
preObj <- preProcess(training[,-58],method=c("center","scale"))
# normalize training set
trainCapAveS <- predict(preObj,training[,-58])$capitalAve
# normalize test set using training parameters
testCapAveS <- predict(preObj,testing[,-58])$capitalAve
# compare results for capitalAve variable
rbind(train = c(mean = mean(trainCapAveS), std = sd(trainCapAveS)),
    test = c(mean(testCapAveS), sd(testCapAveS)))

##               mean      std
## train 6.362510e-18 1.000000
## test  7.548133e-02 1.633866

• preprocess(data, method="BoxCox") = applies BoxCox transformations to continuous data to help normalize the variables through maximum likelihood
  • Note: note it assumes continuous values and DOES NOT deal with repeated values
  • qqnorm(processedVar) = can be used to produce the Q-Q plot which compares the theoretical quantiles with the sample quantiles to see the normality of the data

library(e1071)
# set up BoxCox transforms
preObj <- preProcess(training[,-58],method=c("BoxCox"))
# perform preprocessing on training data
trainCapAveS <- predict(preObj,training[,-58])$capitalAve
# plot histogram and QQ Plot
# Note: the transformation definitely helped to
# normalize the data but it does not produce perfect result
par(mfrow=c(1,2)); hist(trainCapAveS); qqnorm(trainCapAveS)

• preProcess(data, method="knnImpute") = impute/estimate the missing data using k nearest neighbors (knn) imputation
  • knnImpute = takes the k nearest neighbors from the missing value and averages the value to impute the missing observations
  • Note: most prediction algorithms are not build to handle missing data

# Make some values NA
training$capAve <- training$capitalAve
selectNA <- rbinom(dim(training)[1],size=1,prob=0.05)==1
training$capAve[selectNA] <- NA
# Impute and standardize
preObj <- preProcess(training[,-58],method="knnImpute")
capAve <- predict(preObj,training[,-58])$capAve
# Standardize true values
capAveTruth <- training$capitalAve
capAveTruth <- (capAveTruth-mean(capAveTruth))/sd(capAveTruth)
# compute differences between imputed values and true values
quantile(capAve - capAveTruth)

##            0%           25%           50%           75%          100% 
## -1.8242434686 -0.0005069839  0.0007532392  0.0013598891  0.4062295931

Covariate Creation/Feature Extraction

• [level 1]: construct covariate (usable metric, feature) from raw data depends heavily on application
  • ideally we want to summarize data without too much information loss
  • examples
    • text files: frequency of words, frequency of phrases (Google ngrams), frequency of capital letters
    • images: edges, corners, blobs, ridges (computer vision feature detection)
    • webpages: number and type of images, position of elements, colors, videos (A/B Testing)
    • people: height, weight, hair color, sex, country of origin
  • generally, more knowledge and understanding you have of the system/data, the easier it will be to extract the summarizing features
    • when in doubt, more features is always safer –> lose less information and the features can be filtered during model construction
  • this process can be automated (i.e. PCA) but generally have to be very careful, as one very useful feature in the training data set may not have as much effect on the test data set
  • Note: science is the key here, Google “feature extraction for [data type]” for more guidance
    • the goal is always to find the salient characteristics that are likely to be different from observation to observation
• [level 2]: construct new covariates from extracted covariate
  • generally transformations of features you extract from raw data
  • used more for methods like regression and support vector machines (SVM), whose accuracy depend more on the distribution of input variables
  • models like classification trees don’t require as many complex covariates
  • best approach is through exploratory analysis (tables/plots)
  • should only be performed on the train dataset
  • new covariates should be added to data frames under recognizable names so they can be used later
  • preProcess() can be leveraged to handle creating new covariates
  • Note: always be careful about over-fitting

Creating Dummy Variables

• convert factor variables to indicator/dummy variable –> qualitative become quantitative
• dummyVars(outcome~var, data=training) = creates a dummy variable object that can be used through predict function to create dummy variables
  • predict(dummyObj, newdata=training) = creates appropriate columns to represent the factor variable with appropriate 0s and 1s
    • 2 factor variable –> two columns which have 0 or 1 depending on the outcome
    • 3 factor variable –> three columns which have 0, 0, and 1 representing the outcome
    • Note: only one of the columns can have values of 1 for each observation

# setting up data
inTrain <- createDataPartition(y=Wage$wage,p=0.7, list=FALSE)
training <- Wage[inTrain,]; testing <- Wage[-inTrain,]
# create a dummy variable object
dummies <- dummyVars(wage ~ jobclass,data=training)
# create the dummy variable columns
head(predict(dummies,newdata=training))

##        jobclass.1. Industrial jobclass.2. Information
## 231655                      1                       0
## 86582                       0                       1
## 11443                       0                       1
## 305706                      1                       0
## 160191                      1                       0
## 448410                      1                       0

Removing Zero Covariates

• some variables have no variability at all (i.e. variable indicating if an email contained letters)
• these variables are not useful when we want to construct a prediction model
• nearZeroVar(training, saveMetrics=TRUE) = returns list of variables in training data set with information on frequency ratios, percent uniques, whether or not it has zero variance
  • freqRatio = ratio of frequencies for the most common value over second most common value
  • percentUnique = percentage of unique data points out of total number of data points
  • zeroVar = TRUE/FALSE indicating whether the predictor has only one distinct value
  • nzv = TRUE/FALSE indicating whether the predictor is a near zero variance predictor
  • Note: when nzv = TRUE, those variables should be thrown out

# print nearZeroVar table
nearZeroVar(training,saveMetrics=TRUE)

##            freqRatio percentUnique zeroVar   nzv
## year        1.017647    0.33301618   FALSE FALSE
## age         1.231884    2.85442436   FALSE FALSE
## sex         0.000000    0.04757374    TRUE  TRUE
## maritl      3.329571    0.23786870   FALSE FALSE
## race        8.480583    0.19029496   FALSE FALSE
## education   1.393750    0.23786870   FALSE FALSE
## region      0.000000    0.04757374    TRUE  TRUE
## jobclass    1.070936    0.09514748   FALSE FALSE
## health      2.526846    0.09514748   FALSE FALSE
## health_ins  2.209160    0.09514748   FALSE FALSE
## logwage     1.011765   18.83920076   FALSE FALSE
## wage        1.011765   18.83920076   FALSE FALSE

Creating Splines (Polynomial Functions)

• when you want to fit curves through the data, basis functions can be leveraged
• [splines package] bs(data$var, df=3) = creates 3 new columns corresponding to the var, var², and var³ terms
  
• ns() and poly() can also be used to generate polynomials
• gam() function can also be used and it allows for smoothing of multiple variables with different values for each variable
• Note: the same polynomial operations must be performed on the test sets using the predict function

# load splines package
library(splines)
# create polynomial function
bsBasis <- bs(training$age,df=3)
# fit the outcome on the three polynomial terms
lm1 <- lm(wage ~ bsBasis,data=training)
# plot all age vs wage data
plot(training$age,training$wage,pch=19,cex=0.5)
# plot the fitted polynomial function
points(training$age,predict(lm1,newdata=training),col="red",pch=19,cex=0.5)

# predict on test values
head(predict(bsBasis,age=testing$age))

##               1          2           3
## [1,] 0.00000000 0.00000000 0.000000000
## [2,] 0.23685006 0.02537679 0.000906314
## [3,] 0.36252559 0.38669397 0.137491189
## [4,] 0.42616898 0.14823269 0.017186399
## [5,] 0.44221829 0.19539878 0.028779665
## [6,] 0.01793746 0.20448709 0.777050955

Multicore Parallel Processing

• many of the algorithms in the caret package are computationally intensive
• since most of the modern machines have multiple cores on their CPUs, it is often wise to enable multicore parallel processing to expedite the computations
• doMC package is recommended to be used for caret computations (reference)
  • doMC::registerDoMC(cores=4) = registers 4 cores for R to utilize
  • the number of cores you should specify depends on the CPU on your computer (system information usually contains the number of cores)
    • it’s also possible to find the number of cores by directly searching for your CPU model number on the Internet
  • Note: once registered, you should see in your task manager/activity monitor that 4 “R Session” appear when you run your code

Preprocessing with Principal Component Analysis (PCA)

• constructing a prediction model may not require every predictor
• ideally we want to capture the most variation with the least amount of variables
  • weighted combination of predictors may improve fit
  • combination needs to capture the most information
• PCA is suited to do this and will help reduce number of predictors as well as reduce noise (due to averaging)
  • statistical goal = find new set of multivariate variables that are uncorrelated and explain as much variance as possible
  • data compression goal = find the best matrix created with fewer variables that explains the original data
  • PCA is most useful for linear-type models (GLM, LDA)
  • generally more difficult to interpret the predictors (complex weighted sums of variables)
  • Note: outliers are can be detrimental to PCA as they may represent a lot of variation in data
    • exploratory analysis (plots/tables) should be used to identify problems with the predictors
    • transformations with log/BoxCox may be helpful

prcomp Function

• pr<-prcomp(data) = performs PCA on all variables and returns a prcomp object that contains information about standard deviations and rotations
  • pr$rotations = returns eigenvectors for the linear combinations of all variables (coefficients that variables are multiplied by to come up with the principal components) –> how the principal components are created
  • often times, it is useful to take the log transformation of the variables and adding 1 before performing PCA
    • helps to reduce skewness or strange distribution in data
    • log(0) = - infinity, so we add 1 to account for zero values
    • makes data more Gaussian
  • plot(pr) = plots the percent variation explained by the first 10 principal components (PC)
    • can be used to find the PCs that represent the most variation

# load  spam data
data(spam)
# perform PCA on dataset
prComp <- prcomp(log10(spam[,-58]+1))
# print out the eigenvector/rotations first 5 rows and PCs
head(prComp$rotation[, 1:5], 5)

##                 PC1           PC2         PC3         PC4          PC5
## make    0.019370409  0.0427855959 -0.01631961  0.02798232 -0.014903314
## address 0.010827343  0.0408943785  0.07074906 -0.01407049  0.037237531
## all     0.040923168  0.0825569578 -0.03603222  0.04563653  0.001222215
## num3d   0.006486834 -0.0001333549  0.01234374 -0.01005991 -0.001282330
## our     0.036963221  0.0941456085 -0.01871090  0.05098463 -0.010582039

# create new variable that marks spam as 2 and nospam as 1
typeColor <- ((spam$type=="spam")*1 + 1)
# plot the first two principal components
plot(prComp$x[,1],prComp$x[,2],col=typeColor,xlab="PC1",ylab="PC2")

caret Package

• pp<-preProcess(log10(training[,-58]+1),method="pca",pcaComp=2,thresh=0.8)) = perform PCA with preProcess function and returns the number of principal components that can capture the majority of the variation
  • creates a preProcess object that can be applied using predict function
  • pcaComp=2 = specifies the number of principal components to compute (2 in this case)
  • thresh=0.8 = threshold for variation captured by principal components
    • thresh=0.95 = default value, which returns the number of principal components that are needed to capture 95% of the variation in data
• predict(pp, training) = computes new variables for the PCs (2 in this case) for the training data set
  • the results from predict can then be used as data for the prediction model
  • Note: the same PCA must be performed on the test set

# create train and test sets
inTrain <- createDataPartition(y=spam$type,p=0.75, list=FALSE)
training <- spam[inTrain,]
testing <- spam[-inTrain,]
# create preprocess object
preProc <- preProcess(log10(training[,-58]+1),method="pca",pcaComp=2)
# calculate PCs for training data
trainPC <- predict(preProc,log10(training[,-58]+1))
# run model on outcome and principle components
modelFit <- train(training$type ~ .,method="glm",data=trainPC)
# calculate PCs for test data
testPC <- predict(preProc,log10(testing[,-58]+1))
# compare results
confusionMatrix(testing$type,predict(modelFit,testPC))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction nonspam spam
##    nonspam     656   41
##    spam         82  371
##                                           
##                Accuracy : 0.893           
##                  95% CI : (0.8737, 0.9103)
##     No Information Rate : 0.6417          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.7724          
##  Mcnemar's Test P-Value : 0.0003101       
##                                           
##             Sensitivity : 0.8889          
##             Specificity : 0.9005          
##          Pos Pred Value : 0.9412          
##          Neg Pred Value : 0.8190          
##              Prevalence : 0.6417          
##          Detection Rate : 0.5704          
##    Detection Prevalence : 0.6061          
##       Balanced Accuracy : 0.8947          
##                                           
##        'Positive' Class : nonspam         
## 

• alternatively, PCA can be directly performed with the train method
  • train(outcome ~ ., method="glm", preProcess="pca", data=training) = performs PCA first on the training set and then runs the specified model
    • effectively the same procedures as above (preProcess –> predict)
• Note: in both cases, the PCs were able to achieve 90+% accuracy

# construct model
modelFit <- train(training$type ~ .,method="glm",preProcess="pca",data=training)
# print results of model
confusionMatrix(testing$type,predict(modelFit,testing))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction nonspam spam
##    nonspam     668   29
##    spam         59  394
##                                           
##                Accuracy : 0.9235          
##                  95% CI : (0.9066, 0.9382)
##     No Information Rate : 0.6322          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.8379          
##  Mcnemar's Test P-Value : 0.001992        
##                                           
##             Sensitivity : 0.9188          
##             Specificity : 0.9314          
##          Pos Pred Value : 0.9584          
##          Neg Pred Value : 0.8698          
##              Prevalence : 0.6322          
##          Detection Rate : 0.5809          
##    Detection Prevalence : 0.6061          
##       Balanced Accuracy : 0.9251          
##                                           
##        'Positive' Class : nonspam         
## 

Predicting with Regression

• prediction with regression = fitting regression model (line) to data –> multiplying each variable by coefficients to predict outcome
• useful when the relationship between the variables can be modeled as linear
• the model is easy to implement and the coefficients are easy to interpret
• if the relationships are non-linear, the regression model may produce poor results/accuracy
  • Note: linear regressions are generally used in combination with other models
• model \[Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \ldots + \beta_p X_{pi} + e_i \]
  • where \(\beta_0\) is the intercept (when all variables are 0)
  • \(\beta_1, \beta_2, \ldots, \beta_p\) are the coefficients
  • \(X_{1i}, X_{2i}, \ldots, X_{pi}\) are the variables/covariates
  • \(e_i\) is the error
  • \(Y_i\) is the outcome
• prediction \[\hat Y_i = \hat \beta_0 + \hat \beta_1 X_{1i} + \hat \beta_2 X_{2i} + \ldots + \hat \beta_p X_{pi}\]
  • where \(\hat \beta_0\) is the estimated intercept (when all variables are 0)
  • \(\hat \beta_1, \hat \beta_2, \ldots, \hat \beta_p\) are the estimated coefficients
  • \(X_{1i}, X_{2i}, \ldots, X_{pi}\) are the variables/covariates
  • \(\hat Y_i\) is the predicted outcome

R Commands and Examples

• lm<-lm(y ~ x, data=train) = runs a linear model of outcome y on predictor x –> univariate regression
  • summary(lm) = returns summary of the linear regression model, which will include coefficients, standard errors, \(t\) statistics, and p values
  • lm(y ~ x1+x2+x3, data=train) = run linear model of outcome y on predictors x1, x2, and x3
  • lm(y ~ ., data=train = run linear model of outcome y on all predictors
• predict(lm, newdata=df) = use the constructed linear model to predict outcomes (\(\hat Y_i\))for the new values
  • newdata data frame must have the same variables (factors must have the same levels) as the training data
  • newdata=test = predict outcomes for the test set based on linear regression model from the training
  • Note: the regression line will not be a perfect fit on the test set since it was constructed on the training set
• RSME can be calculated to measure the accuracy of the linear model
  • Note: \(RSME_{test}\), which estimates the out-of-sample error, is almost always GREATER than \(RSME_{train}\)

# load data
data(faithful)
# create train and test sets
inTrain <- createDataPartition(y=faithful$waiting, p=0.5, list=FALSE)
trainFaith <- faithful[inTrain,]; testFaith <- faithful[-inTrain,]
# build linear model
lm1 <- lm(eruptions ~ waiting,data=trainFaith)
# print summary of linear model
summary(lm1)

## 
## Call:
## lm(formula = eruptions ~ waiting, data = trainFaith)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.30541 -0.35970  0.04711  0.37624  1.07683 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.870699   0.225663   -8.29 1.01e-13 ***
## waiting      0.075674   0.003137   24.12  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5024 on 135 degrees of freedom
## Multiple R-squared:  0.8117, Adjusted R-squared:  0.8103 
## F-statistic: 581.8 on 1 and 135 DF,  p-value: < 2.2e-16

# predict eruptions for new waiting time
newdata <- data.frame(waiting=80)
predict(lm1,newdata)

##        1 
## 4.183192

# create 1 x 2 panel plot
par(mfrow=c(1,2))
# plot train data with the regression line
plot(trainFaith$waiting,trainFaith$eruptions,pch=19,col="blue",xlab="Waiting",
    ylab="Duration", main = "Train")
lines(trainFaith$waiting,predict(lm1),lwd=3)
# plot test data with the regression line
plot(testFaith$waiting,testFaith$eruptions,pch=19,col="blue",xlab="Waiting",
    ylab="Duration", main = "Test")
lines(testFaith$waiting,predict(lm1,newdata=testFaith),lwd=3)

# Calculate RMSE on training and test sets
c(trainRMSE = sqrt(sum((lm1$fitted-trainFaith$eruptions)^2)),
    testRMSE = sqrt(sum((predict(lm1,newdata=testFaith)-testFaith$eruptions)^2)))

## trainRMSE  testRMSE 
##  5.836980  5.701161

• pi<-predict(lm, newdata=test, interval="prediction") = returns 3 columns for fit (predicted value, same as before), lwr (lower bound of prediction interval), and upr (upper bound of prediction interval)
  • matlines(x, pi, type="l") = plots three lines, one for the linear fit and two for upper/lower prediction interval bounds

# calculate prediction interval
pred1 <- predict(lm1,newdata=testFaith,interval="prediction")
# plot data points (eruptions, waiting)
plot(testFaith$waiting,testFaith$eruptions,pch=19,col="blue")
# plot fit line and prediction interval
matlines(testFaith$waiting,pred1,type="l",,col=c(1,2,2),lty = c(1,1,1), lwd=3)

• lm <- train(y ~ x, method="lm", data=train) = run linear model on the training data –> identical to lm function
  • summary(lm$finalModel) = returns summary of the linear regression model, which will include coefficients, standard errors, \(t\) statistics, and p values –> identical to summary(lm) for a lm object
  • train(y ~ ., method="lm", data=train) = run linear model on all predictors in training data
    • multiple predictors (dummy/indicator variables) are created for factor variables
  • plot(lm$finalModel) = construct 4 diagnostic plots for evaluating the model
    • Note: more information on these plots can be found at ?plot.lm
    • Residual vs Fitted
    • Normal Q-Q
    • Scale-Location
    • Residual vs Leverage

# create train and test sets
inTrain <- createDataPartition(y=Wage$wage,p=0.7, list=FALSE)
training <- Wage[inTrain,]; testing <- Wage[-inTrain,]
# fit linear model for age jobclass and education
modFit<- train(wage ~ age + jobclass + education,method = "lm",data=training)
# store final model
finMod <- modFit$finalModel
# set up 2 x 2 panel plot
par(mfrow = c(2, 2))
# construct diagnostic plots for model
plot(finMod,pch=19,cex=0.5,col="#00000010")

• plotting residuals by index can be helpful in showing missing variables
  • plot(finMod$residuals) = plot the residuals against index (row number)
  • if there’s a trend/pattern in the residuals, it is highly likely that another variable (such as age/time) should be included
    • residuals should not have relationship to index

# plot residual by index
plot(finMod$residuals,pch=19,cex=0.5)

* here the residuals increase linearly with the index, and the highest residuals are concentrated in the higher indexes, so there must be a missing variable

Prediction with Trees

• prediction with trees = iteratively split variables into groups (effectively constructing decision trees) –> produces nonlinear model
  • the classification tree uses interactions between variables –> the ultimate groups/leafs may depend on many variables
• the result (tree) is easy to interpret, and generally performs better predictions than regression models when the relationships are non-linear
• transformations less important –> monotone transformations (order unchanged, such as \(\log\)) will produce same splits
• trees can be used for regression problems as well and use RMSE as measure of impurity
• however, without proper cross-validation, the model can be over-fitted (especially with large number of variables) and results may be variable from one run to the next
  • it is also harder to estimate the uncertainty of the model
• party, rpart, tree packages can all build trees

Process

1. start with all variables in one group
2. find the variable that best splits the outcomes into two groups
3. divide data into two groups (leaves) based on the split performed (node)
4. within each split, find variables to split the groups again
5. continue this process until all groups are sufficiently small/homogeneous/“pure”

Measures of Impurity (Reference)

\[\hat{p}_{mk} = \frac{\sum_{i}^m \mathbb{1}(y_i = k)}{N_m}\]

• \(\hat p_mk\) is the probability of the objects in group \(m\) to take on the classification \(k\)

• \(N_m\) is the size of the group

• Misclassification Error \[ 1 - \hat{p}_{m~k(m)}\] where \(k(m)\) is the most common classification/group
  • 0 = perfect purity
  • 0.5 = no purity
    • Note: it is not 1 here because when \(\hat{p}_{m~k(m)}\) < 0.5 or there’s not predominant classification for the objects, it means the group should be further subdivided until there’s a majority
• Gini Index\[ \sum_{k \neq k'} \hat{p}_{mk} \times \hat{p}_{mk'} = \sum_{k=1}^K \hat{p}_{mk}(1-\hat{p}_{mk}) = 1 - \sum_{k=1}^K p_{mk}^2\]
  • 0 = perfect purity
  • 0.5 = no purity
• Deviance \[ -\sum_{k=1}^K \hat{p}_{mk} \log_e\hat{p}_{mk} \]
  • 0 = perfect purity
  • 1 = no purity
• Information Gai \[ -\sum_{k=1}^K \hat{p}_{mk} \log_2\hat{p}_{mk} \]
  • 0 = perfect purity
  • 1 = no purity

• example

# set margin and seed
par(mar=c(1,1,1,1), mfrow = c(1, 2)); set.seed(1234);
# simulate data
x = rep(1:4,each=4); y = rep(1:4,4)
# plot first scenario
plot(x,y,xaxt="n",yaxt="n",cex=3,col=c(rep("blue",15),rep("red",1)),pch=19)
# plot second scenario
plot(x,y,xaxt="n",yaxt="n",cex=3,col=c(rep("blue",8),rep("red",8)),pch=19)

* left graph - Misclassification: \(\frac{1}{16} = 0.06\) - Gini: \(1 - [(\frac{1}{16})^2 + (\frac{15}{16})^2] = 0.12\) - Information:\(-[\frac{1}{16} \times \log_2 (\frac{1}{16})+ \frac{15}{16} \times \log_2(\frac{1}{16})] = 0.34\) * right graph - Misclassification: \(\frac{8}{16} = 0.5\) - Gini: \(1 - [(\frac{8}{16})^2 + (\frac{8}{16})^2] = 0.5\) - Information:\(-[\frac{8}{16} \times \log_2 (\frac{8}{16})+ \frac{8}{16} \times \log_2(\frac{8}{16})] = 1\)

Constructing Trees with caret Package

• tree<-train(y ~ ., data=train, method="rpart") = constructs trees based on the outcome and predictors
  • produces an rpart object, which can be used to predict new/test values
  • print(tree$finalModel) = returns text summary of all nodes/splits in the tree constructed
• plot(tree$finalModel, uniform=TRUE) = plots the classification tree with all nodes/splits
  • [rattle package] fancyRpartPlot(tree$finalModel) = produces more readable, better formatted classification tree diagrams
  • each split will have the condition/node in bold and the splits/leafs on the left and right sides following the “yes” or “no” indicators
    • “yes” –> go left
    • “no” –> go right

library(rattle)
library(rpart)
# load iris data set
data(iris)
# create test/train data sets
inTrain <- createDataPartition(y=iris$Species,p=0.7, list=FALSE)
training <- iris[inTrain,]
testing <- iris[-inTrain,]
# fit classification tree as a model
modFit <- train(Species ~ .,method="rpart",data=training)
# print the classification tree
print(modFit$finalModel)

## n= 105 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
## 1) root 105 70 setosa (0.33333333 0.33333333 0.33333333)  
##   2) Petal.Length< 2.45 35  0 setosa (1.00000000 0.00000000 0.00000000) *
##   3) Petal.Length>=2.45 70 35 versicolor (0.00000000 0.50000000 0.50000000)  
##     6) Petal.Width< 1.65 34  1 versicolor (0.00000000 0.97058824 0.02941176) *
##     7) Petal.Width>=1.65 36  2 virginica (0.00000000 0.05555556 0.94444444) *

# plot the classification tree
rattle::fancyRpartPlot(modFit$finalModel)

# predict on test values
predict(modFit,newdata=testing)

##  [1] setosa     setosa     setosa     setosa     setosa     setosa    
##  [7] setosa     setosa     setosa     setosa     setosa     setosa    
## [13] setosa     setosa     setosa     versicolor versicolor versicolor
## [19] versicolor versicolor versicolor versicolor versicolor versicolor
## [25] versicolor versicolor versicolor versicolor versicolor versicolor
## [31] virginica  virginica  virginica  virginica  virginica  virginica 
## [37] versicolor virginica  virginica  versicolor versicolor virginica 
## [43] virginica  virginica  virginica 
## Levels: setosa versicolor virginica

Bagging

• bagging = bootstrap aggregating
  • resample training data set (with replacement) and recalculate predictions
  • average the predictions together or majority vote
  • more information can be found here
• averaging multiple complex models have similar bias as each of the models on its own, and reduced variance because of the average

• most useful for non-linear models

• example
  • loess(y ~ x, data=train, span=0.2) = fits a smooth curve to data
    • span=0.2 = controls how smooth the curve should be

# load data
library(ElemStatLearn); data(ozone,package="ElemStatLearn")
# reorder rows based on ozone variable
ozone <- ozone[order(ozone$ozone),]
# create empty matrix
ll <- matrix(NA,nrow=10,ncol=155)
# iterate 10 times
for(i in 1:10){
    # create sample from data with replacement
    ss <- sample(1:dim(ozone)[1],replace=T)
    # draw sample from the dataa and reorder rows based on ozone
    ozone0 <- ozone[ss,]; ozone0 <- ozone0[order(ozone0$ozone),]
    # fit loess function through data (similar to spline)
    loess0 <- loess(temperature ~ ozone,data=ozone0,span=0.2)
    # prediction from loess curve for the same values each time
    ll[i,] <- predict(loess0,newdata=data.frame(ozone=1:155))
}
# plot the data points
plot(ozone$ozone,ozone$temperature,pch=19,cex=0.5)
# plot each prediction model
for(i in 1:10){lines(1:155,ll[i,],col="grey",lwd=2)}
# plot the average in red
lines(1:155,apply(ll,2,mean),col="red",lwd=2)

Bagging Algorithms

• in the caret package, there are three options for the train function to perform bagging
  • bagEarth - Bagged MARS (documentation)
  • treebag - Bagged CART (documentation)
  • bagFDA - Bagged Flexible Discriminant Analysis (documentation)
• alternatively, custom bag functions can be constructed (documentation)
  • bag(predictors, outcome, B=10, bagControl(fit, predict, aggregate)) = define and execute custom bagging algorithm
    • B=10 = iterations/resampling to perform
    • bagControl() = controls for how the bagging should be executed
      • fit=ctreeBag$fit = the model ran on each resampling of data
      • predict=ctreeBag$predict = how predictions should be calculated from each model
      • aggregate=ctreeBag$aggregate = how the prediction models should be combined/averaged
• example

# load relevant package and data
library(party); data(ozone,package="ElemStatLearn")
# reorder rows based on ozone variable
ozone <- ozone[order(ozone$ozone),]
# extract predictors
predictors <- data.frame(ozone=ozone$ozone)
# extract outcome
temperature <- ozone$temperature
# run bagging algorithm
treebag <- bag(predictors, temperature, B = 10,
                # custom bagging function
                bagControl = bagControl(fit = ctreeBag$fit,
                                        predict = ctreeBag$pred,
                                        aggregate = ctreeBag$aggregate))
# plot data points
plot(ozone$ozone,temperature,col='lightgrey',pch=19)
# plot the first fit
points(ozone$ozone,predict(treebag$fits[[1]]$fit,predictors),pch=19,col="red")
# plot the aggregated predictions
points(ozone$ozone,predict(treebag,predictors),pch=19,col="blue")

Random Forest

• random forest = extension of bagging on classification/regression trees
  • one of the most used/accurate algorithms along with boosting

* process - bootstrap samples from training data (with replacement) - split and bootstrap variables - grow trees (repeat split/bootstrap) and vote/average final trees

• drawbacks
  • algorithm can be slow (process large number of trees)
  • hard to interpret (large numbers of splits and nodes)
  • over-fitting (difficult to know which tree is causing over-fitting)
  • Note: it is extremely important to use cross validation when running random forest algorithms

R Commands and Examples

• rf<-train(outcome ~ ., data=train, method="rf", prox=TRUE, ntree=500) = runs random forest algorithm on the training data against all predictors
  • Note: random forest algorithm automatically bootstrap by default, but it is still important to have train/test/validation split to verify the accuracy of the model
  • prox=TRUE = the proximity measures between observations should be calculated (used in functions such as classCenter() to find center of groups)
    • rf$finalModel$prox = returns matrix of proximities
  • ntree=500 = specify number of trees that should be constructed
  • do.trace=TRUE = prints logs as the trees are being built –> useful by indicating progress to user
  • Note: randomForest() function can be used to perform random forest algorithm (syntax is the same as train) and is much faster
• getTree(rf$finalModel, k=2) = return specific tree from random forest model
• classCenters(predictors, outcome, proximity, nNbr) = return computes the cluster centers using the nNbr nearest neighbors of the observations
  • prox = rf$finalModel$prox = proximity matrix from the random forest model
  • nNbr = number of nearest neighbors that should be used to compute cluster centers
• predict(rf, test) = apply the random forest model to test data set
  • confusionMatrix(predictions, actualOutcome) = tabulates the predictions of the model against the truths
    • Note: this is generally done for the validation data set using the model built from training
• example

library(randomForest)
# load data
data(iris)
# create train/test data sets
inTrain <- createDataPartition(y=iris$Species,p=0.7, list=FALSE)
training <- iris[inTrain,]
testing <- iris[-inTrain,]
# apply random forest
modFit <- train(Species~ .,data=training,method="rf",prox=TRUE)
# return the second tree (first 6 rows)
head(getTree(modFit$finalModel,k=2))

##   left daughter right daughter split var split point status prediction
## 1             2              3         3        2.60      1          0
## 2             0              0         0        0.00     -1          1
## 3             4              5         4        1.65      1          0
## 4             6              7         3        5.25      1          0
## 5             0              0         0        0.00     -1          3
## 6             0              0         0        0.00     -1          2

# compute cluster centers
irisP <- classCenter(training[,c(3,4)], training$Species, modFit$finalModel$prox)
# convert irisP to data frame and add Species column
irisP <- as.data.frame(irisP); irisP$Species <- rownames(irisP)
# plot data points
p <- qplot(Petal.Width, Petal.Length, col=Species,data=training)
# add the cluster centers
p + geom_point(aes(x=Petal.Width,y=Petal.Length,col=Species),size=5,shape=4,data=irisP)

# predict outcome for test data set using the random forest model
pred <- predict(modFit,testing)
# logic value for whether or not the rf algorithm predicted correctly
testing$predRight <- pred==testing$Species
# tabulate results
table(pred,testing$Species)

##             
## pred         setosa versicolor virginica
##   setosa         15          0         0
##   versicolor      0         13         1
##   virginica       0          2        14

# plot data points with the incorrect classification highlighted
qplot(Petal.Width,Petal.Length,colour=predRight,data=testing,main="newdata Predictions")

Boosting

• boosting = one of the most widely used and accurate prediction models, along with random forest
• boosting can be done with any set of classifiers, and a well-known approach is gradient boosting

• more detail tutorial can be found here

• process: take a group of weak predictors –> weight them and add them up –> result in a stronger predictor
  • start with a set of classifiers \(h_1, \ldots, h_k\)
    • examples: all possible trees, all possible regression models, all possible cutoffs (divide data into different parts)
  • calculate a weighted sum of classifiers as the prediction value \[f(x) = \sum_i \alpha_i h_i(x)\] where \(\alpha_i\) = coefficient/weight and \(h_i(x)\) = value of classifier
    • goal = minimize error (on training set)
    • select one \(h\) at each step (iterative)
    • calculate weights based on errors
    • up-weight missed classifications and select next \(h\)

• example

* we start with space with blue + and red - and the goal is to classify all the object correctly * only straight lines will be used for classification

* from the above, we can see that a group of weak predictors (lines in this case), can be combined and weighed to become a much stronger predictor